% =========================================================================
% lmj_parperturb
% 1. Load in
% a) one  parameter mode and perturb one coefficient at a time 
%    according to partvec. True value will be added automatically 
% b) a table of modes 
% Also load the dataset 
% =========================================================================
clear all;close all;clc;cucd=pwd; osncreen=0; 

%% USER INPUTS
loadpath = 'O:\PROJ_LIB\LJM\exog\udsregmebnd\lindtr\1trend z terragona mode'; 
presetf = 'preset_eval'; 
compare_type = 2; 
savefolder='terragona '; 
par_pert = 'alpha';
pertvec = [0.15 0.50];
modenames=[]; 

%% 
% Load in settings file
cd(loadpath); 
eval(presetf); 
cd(cucd); 

% Define Save_path
if isempty(savefolder) && compare_type == 1; 
   savefolder = ['par_perturb', strdate, '\', par_pert]; 
elseif compare_type == 1
    savefolder = [savefolder, strdate, '\', par_pert]; 
else 
   savefolder = [savefolder, strdate]; 
end
save_path = cr_dir(loadpath,savefolder); 


% Load in Parameter Modes
cd(loadpath); 
[parmode, parnames] = xlsread(parf_name, partab_name); 
cd(cucd); 
numpar_loaded = size(parmode, 2);


% Define function handle
temp      = extrsubf(loadpath, cucd ); 
tempos    = findstr(temp,'\'); 
root=temp(tempos(1):tempos(2)-1); 
spec=temp(tempos(2):tempos(3)-1); 
subf=temp(tempos(3):end        ); 
clear temp tempos; 
funcmod=['mod',root(2:end), spec(2:end)]; 

% =========================================================================
% Change function handle to FUNCMOD + ADD2FUNC 
% =========================================================================
if exist('add2func','var')==1 && ~isempty(add2func) 
    dispaj('Switching function handle from ',funcmod)
    dispaj('to  ',[funcmod,add2func]);
    disp(' ');
    disp('Verifying solutions match...');
    pause(0.25);
    [var1,var2,difvar,difc]=checkmodels_lmj(str2func(funcmod),...
        str2func([funcmod,add2func]),parmode(:,1),solveopt,addsol);
    if difvar < 1e-7;
        disp('Solutions identical')
    else
        disp('Solutions do not match')
    end
    
    if difc < 1e-10;
        disp('Constants identical');
    else
        disp('Constants do not match');
    end
    clear var1 var2 difvar difc
    funcmod=str2func([funcmod,add2func]);
    
    disp('Changed funcmod handle')        
else 
    funcmod=str2func(funcmod); 
end 


% check things
if compare_type == 1 
    if numpar_loaded ~= 1; 
        error('only one parameter vector allowed under PERTURB mode');
    elseif ~exist('par_pert', 'var') || isempty('par_pert');  
        error('perturbed variable not selected under PERTURB mode');
    elseif ~exist('pertvec', 'var') || isempty('pertvec');
        error('PERTVEC not defined under PERTURB mode'); 
    end
end 

% if compare_type == 2 && numpar_loaded == 1; 
%     error('at leaset two parameter vector required under COMPARE mode'); 
% end 


printcell([parnames num2cprec(parmode)]); 
if onscreen==1
    quer('c');
end

% Load in Data
cd(dataf_path); 
[Y,Ynames,sample,trainvec]=load_dset(dataf_name,dataf_path,dataopt,dataopt.demean_type,onscreen);
cd(cucd); 

% check each loaded modes yields good results
for jj = 1:numpar_loaded
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
        =feval(funcmod,parmode(:, jj),solveopt,addsol);
    if ~isequal(eu, [1; 1]);
        error('loaded parameter does not solve the model ')
    end
    znames=cell(size(ZZ,1),1);
    for ii=1:length(znames)
        temp      =find(ZZ(ii,:)==1 ) ;
        if length(temp) > 1
            error('Need 1 entry only for each row of ZZ')
        end
        znames(ii)=stanames( temp );
    end
end


%% solve for the model
%  Sorted in ascending order 

addmodenames = isempty(modenames); 

switch compare_type
    case 1
        parpos_pert = cellposition(par_pert, parnames);
        parpert_true = parmode(parpos_pert);
        pos_dupvar = find(abs(parpert_true - pertvec) < 1e-5);
        if ~isempty(pos_dupvar)
            pertvec(pos_dupvar) = nan;
            pertvec = pertvec(~isnan(pertvec));
        end
        pertvec=[parpert_true, pertvec];
        numpert = length(pertvec);
        pertmat = repmat(parmode, 1, numpert);
        pertmat(parpos_pert,:) = pertvec;
        % Create top cell
        modenames = cell(1, numpert);
        for ii = 1:numpert;
            if abs( pertvec(ii) - parpert_true ) < eps
                modenames{ii} = [par_pert, ' = ', num2str(parpert_true), ' (est)'];
            else
                modenames{ii} = [par_pert, ' = ', num2str(pertvec(ii))];
            end
        end
    case 2        
        % Check if modenames exists, assing default otherwise 
        pertmat = parmode;
        numpert = numpar_loaded;
        
        if exist('modenames','var')==1 && ~isempty(modenames)
            if length(modenames)~=numpert
                
                printcell(modenames)
                error('Modenames does not match number of modes loaded')
            end
        else 
            disp('Assing default names to different modes') 
            modenames = cell(1, numpert);
            for ii = 1:numpert;
                modenames{ii} = ['mode ', num2str(ii)];
            end
        end
        
        if onscreen==1 
            disp('Names of modes') 
            printcell(modenames) 
            quer('c') 
        end 
        
end


ssvec_mat = zeros(length(ssnames), numpar_loaded);

for ii = 1:numpar_loaded
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
        =feval(funcmod,parmode(:, ii),solveopt,addsol);
    
    if ii==1                
        [junk,ana_in.flag_adds]=extractlevels(ana_in.state_sel,stanames); 
    end 
    ssvec_mat(:, ii) = ssvec(:);
    figname = ['mode', num2str(ii), '_specdecom'];
    [out,spec_outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,ana_in.specrep,stanames,shonames,znames,...
        ana_in.lowper,ana_in.highper,ana_in.npoints,save_path, figname);
    spec_outcell = merge_cells({['Mode ', num2str(ii)]}, spec_outcell);
    
    [vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
    var_outcell = crtcell( vdcom, znames , shonames ,['Mode ', num2str(ii),': Stationary Variance']);
    
    if ii == 1;
        spec_outcell_full = spec_outcell;
        var_outcell_full = var_outcell;
    else
        spec_outcell_full = merge_cells( spec_outcell_full, spec_outcell, 2);
        var_outcell_full  = merge_cells( var_outcell_full,   var_outcell, 2);
    end
end

outfolder = extrsubf(save_path, loadpath);
topcell = {'root', root;
    'spec', spec;
    'subf', subf;
    'outfolder', outfolder
    'data', dataf_name};

spec_outcell_full = merge_cells(topcell, spec_outcell_full, 2);
var_outcell_full = merge_cells(topcell, var_outcell_full, 2);

par_cell = [parnames, num2cprec(parmode, 10)];
ss_cell = [ssnames, num2cprec(ssvec_mat)];
mode_cell = merge_cells(par_cell, ss_cell, 2);
mode_cell = merge_cells(topcell, mode_cell, 2);
cd(save_path)
xlswrite('compare',spec_outcell_full,'spectrum');
xlswrite('compare', var_outcell_full,'stationary');
xlswrite('compare', mode_cell,'modes');
cd(cucd)

close all



%%

[e, ana_in] = ch_field(ana_in, 'flag_addz',  zeros(1,size(Y, 2)));
if length(znames) ~= length(ana_in.flag_addz);
    error('ana_in.flag_adds not well defined, check setting file')
end
ana_in.znames = znames; 
ana_in.stanames = stanames;
ana_in.shonames = shonames;
ana_in.parnames = parnames;
ana_in.vecsim = [vecsim_d1, size(Y,1)];
ana_in.savefolder = savefolder;
ana_in.modelnames = modenames;
ana_in.ssnames = ssnames;
ana_in.savepath = save_path;

%% mode_analysis COMPUTE
ana_out = lmj_modeanalysis_compute(funcmod, pertmat, Y, ana_in, solveopt, addsol);
pertmat = ana_out.pertmat;

for ii = 1:ana_out.num_valid;
    stdsimz_ii = squeeze(ana_out.stdsimz_mat(:, :, ii)); 
    std_outcell = playmod_outdraws(stdsimz_ii', znames); 
    std_outcell(2, 1:2) = {'model: ', ana_out.modelnames{ii}}; 
    if ii == 1;
        std_outcell_full = std_outcell; 
    else
        std_outcell_full = merge_cells(std_outcell_full, std_outcell, 2 ); 
    end
end 
std_outcell_full = merge_cells(topcell, std_outcell_full); 

cd(save_path)
xlswrite('compare',std_outcell_full,'std_sim');
cd(cucd)


% 1) plot volatilities
in_stdhist.modenames = ana_out.modelnames;
in_stdhist.varnames  = znames;
in_stdhist.savepath  = ana_in.savepath;
in_stdhist.savename  = 'Z STD';
plot_stdhist(ana_out.stdsimz_mat, Y, in_stdhist)

% 2) plot irfs

% (2.1) plot all observables
in_irf.modenames = ana_out.modelnames;
in_irf.varnames = znames;
in_irf.shonames = shonames;
in_irf.savepath = save_path;
in_irf.savename = 'Z IRF';
in_irf.addlev = ana_in.flag_addz;
[e, in_irf] = ch_field(in_irf, 'flag_irfbysho', 1);
plotmultiple_irfs(ana_out.irfz_mat, in_irf);


% (2.2) plotting selected states.
in_irf.modenames = ana_out.modelnames;
in_irf.varnames =  ana_in.state_sel;
in_irf.shonames = shonames;
in_irf.savepath = save_path;
in_irf.savename = 'Sel States IRF';
in_irf.addlev = ana_in.flag_adds;
[e, in_irf] = ch_field(in_irf, 'flag_irfbysho', 1);

plotmultiple_irfs_split(ana_out.irfs_mat, in_irf);


% 3) plot Variance Decompositions.
% 3.1) of the observables
in_vdec.modenames = ana_out.modelnames;
in_vdec.varnames =  znames;
in_vdec.shonames = shonames;
in_vdec.savepath = save_path;
in_vdec.savename = 'Z VDEC';
plot_vdecs(ana_out.vdechz_mat, in_vdec);


% 3.2) of the variables
in_vdec.modenames = ana_out.modelnames;
in_vdec.varnames =  ana_in.state_sel;
in_vdec.shonames = shonames;
in_vdec.savepath = save_path;
in_vdec.savename = 'Sel State VDEC';
plot_vdecs(ana_out.vdechs_mat, in_vdec);


% 4) AUTOCORRELATION
% 4.1) of the obervables;
in_ac.modenames = ['data';  ana_out.modelnames(:)];
in_ac.varnames = znames;
in_ac.savepath = save_path;
in_ac.savename = 'Obs AC Median';

ac_data=accov(Y,(0:ana_in.nper));
clear med_acsimz_mat;
for ii = 1:ana_out.num_valid + 1;
    if ii == 1;
        med_acsimz_mat(:,:,:,ii)=ac_data;
    else
        med_acsimz_mat(:,:,:,ii)=median(ana_out.acsimz_mat(:, :, :, :, ii-1),4);
    end
end
plotmultiple_acorr(med_acsimz_mat, in_ac);

% 4.2) of the selected states
in_ac.modenames = ana_out.modelnames;
in_ac.varnames  = ana_in.state_sel;
in_ac.savepath = save_path;
in_ac.savename = 'Sel States AC Median';

clear med_acsims_mat
for ii = 1:ana_out.num_valid;
    med_acsims_mat(:,:,:,ii)=median(ana_out.acsims_mat(:, :, :, :, ii),4);
end
plotmultiple_acorr(med_acsims_mat, in_ac);


% 4.3) full cross-correlogram vs. data

in_acfull.modenames =  [ana_out.modelnames(:); 'data'];
in_acfull.varnames  = znames;
in_acfull.savepath  = save_path;
in_acfull.savename = 'Obs AC full';
plotmultiple_acall(ana_out.acsimz_mat, ac_data, in_acfull);


% 5, selected steady state values (if perturbing one parameter)
if compare_type == 1
    in_ssplot.modenames = ana_out.modelnames;
    in_ssplot.savepath  = save_path;
    in_ssplot.savename  = 'Sel Steady States';
    in_ssplot.ssreport  = ana_in.ssreport;
    in_ssplot.ssnames   = ssnames;
    in_ssplot.parnames  = parnames;
    in_ssplot.pertvar = par_pert;
    
    plotperturb_lmj(pertmat, funcmod, in_ssplot, solveopt, addsol);
end
